NASA/TP-2012-217573 



Analytic and Computational Perspectives 
of Multi-Scale Theory for Homogeneous, 
Laminated Composite, and Sandwich Beams 
and Plates 


Alexander Tessler 

Langley Research Center, Hampton, Virginia 

Marco Gherlone, Daniele Versino, and Marco Di Sciuva 
Department of Mechanical and Aeronautics Engineering 
Politecnico di Torino, Torino, Italy 


May 2012 




NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA scientific and technical information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI program operates under the 
auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 
disseminates NASA’s STI. The NASA STI 
program provides access to the NASA Aeronautics 
and Space Database and its public interface, the 
NASA Technical Report Server, thus providing one 
of the largest collections of aeronautical and space 
science STI in the world. Results are published in 
both non-NASA channels and by NASA in the 
NASA STI Report Series, which includes the 
following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
Programs and include extensive data or 
theoretical analysis. Includes compilations of 
significant scientific and technical data and 
information deemed to be of continuing 
reference value. NASA counterpart of peer- 
reviewed formal professional papers, but 
having less stringent limitations on manuscript 
length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or of 
specialized interest, e.g., quick release reports, 
working papers, and bibliographies that contain 
minimal annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. 

Collected papers from scientific and 
technical conferences, symposia, seminars, 
or other meetings sponsored or co- 
sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, 
often concerned with subjects having 
substantial public interest. 

• TECHNICAL TRANSLATION. 
English-language translations of foreign 
scientific and technical material pertinent to 
NASA’s mission. 

Specialized services also include organizing 
and publishing research results, distributing 
specialized research announcements and feeds, 
providing information desk and personal search 
support, and enabling data exchange services. 

For more information about the NASA STI 
program, see the following: 

• Access the NASA STI program home page 
at http://www.sti.nasa.sov 

• E-mail your question to help@sti.nasa.gov 

• Fax your question to the NASA STI 
Information Desk at 443-757-5803 

• Phone the NASA STI Information Desk at 
443-757-5802 

• Write to: 

STI Information Desk 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 


NASA/TP-2012-217573 



Analytic and Computational Perspectives 
of Multi-Scale Theory for Homogeneous, 
Laminated Composite, and Sandwich Beams 
and Plates 


Alexander Tessler 

Langley Research Center, Hampton, Virginia 

Marco Gherlone, Daniele Versino, and Marco Di Sciuva 
Department of Mechanical and Aeronautics Engineering 
Politecnico di Torino, Torino, Italy 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virginia 23681-2199 


May 2012 




Available from: 

NASA Center for AeroSpace Information 
7115 Standard Drive 
Hanover, MD 21076-1320 
443-757-5802 





Abstract 


This paper reviews the theoretical foundation and computational 
mechanics aspects of the recently developed shear- deformation theory, 
called the Refined Zigzag Theory (RZT). The theory is based on a multi- 
scale formalism in which an equivalent single-layer plate theory is 
refined with a robust set of zigzag local layer displacements that are free 
of the usual deficiencies found in common plate theories with zigzag 
kinematics. In the RZT. first-order shear-deformation plate theory is 
used as the equivalent single-layer plate theory, which represents the 
overall response characteristics. Local piecewise -linear zigzag 
displacements are used to provide corrections to these overall response 
characteristics that are associated with the plate heterogeneity and the 
relative stiffnesses of the layers. The theory does not rely on shear- 
correction factors and is equally accurate for homogeneous, laminated 
composite, and sandwich beams and plates. Regardless of the number of 
material layers, the theory maintains only seven kinematic unknowns that 
describe the membrane, bending, and transverse shear plate-deformation 
modes. Derived from the virtual work principle, RZT is well-suited for- 
developing computationally efficient, C ° -continuous finite elements; 
formulations of several RZT -based elements are highlighted. The theory 
and its finite element approximations thus provide a unified and reliable 
computational platform for the analysis and design of high-performance 
load-bearing aerospace structures. 
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1 INTRODUCTION 

Lightweight and high-performance characteristics of advanced composite materials have 
spurred a wide range of application of these materials in military and civilian aircraft, aerospace 
vehicles, and naval and civil structures. To reali z e the full potential of composite structures for 
primary load-bearing components, further advances in structural design, analysis methods, and 
failure prediction and progressive damage methodologies are necessary. 

A wide variety of modern civilian and military aircraft use relatively thick laminated 
composite and sandwich laminates for primary load-bearing structures. Such structures 
commonly undergo relatively pronounced design-critical transverse-shear deformations and, 
under certain conditions, thickness-stretch deformations. Fail-safe designs require accurate 
stress-analysis methods, particularly in the regions of stress concentration. Moreover, 
computationally efficient progressive damage analysis necessitates accurate modeling of the 
interlaminar damage modes such as delamination. For these reasons, structural theories that 
account for higher-order deformation effects have attracted much attention in recent years. 

Recently, Tessler and co-authors presented an improved structural theory for beams and 
plates, labeled the Refined Zigzag Theory (RZT), that offers substantial analytic and 
computational advantages for the analysis of homogeneous, laminated composite, and sandwich 
laminates [1-5]. This new theory is based on a multi-scale approach in which an equivalent 
single-layer plate theory is used to represent the overall (coarse) plate response characteristics 
and through-the-thickness zigzag kinematics are used to model the local (fine) layer-level 
behavior. In particular, the RZT uses First-order Shear-Deformation Theory (FSDT) as the 
equivalent single-layer theory, and uses sets of piecewise linear continuous functions to model 
the local behavior of the layers comprising a plate. The zigzag kinematic framework enables 
sufficiently accurate and computationally efficient modeling of a wide range of homogeneous 
and heterogeneous laminates without the use of shear correction factors. Novel zigzag functions, 
derived a priori from constitutive relations without enforcing debilitating stress-equilibrium 
constraints, are responsible for overcoming several critical shortcomings of the earlier zigzag 
theories (refer to [1-5] for the literature reviews and pertinent discussions). The stress resultants 
are obtained from the equilibrium equations and, as a result, are physically consistent with their 
definitions based on Hooke’s relations. The formulation, which maintains a fixed number of 
kinematic unknowns regardless of the number of material layers, does not enforce full continuity 
of the transverse shear stresses along material-layer interfaces, yet is robust. Using the principle 
of virtual work, equilibrium equations and consistent boundary conditions are derived in a 
variationally consistent manner. The variational framework, requiring relatively simple C°- 
continuous kinematic interpolations, provides a convenient means for developing 
computationally efficient and robust finite elements. 

The focus of this paper is to assess the current state of the art in the development of RZT and 
to highlight its recent advances in finite element analysis. To accomplish these objectives, this 
paper first reviews the theoretical foundation of RZT for laminated plates and highlights the 
essential aspects of the piecewise linear zigzag functions derived from transverse-shear 
constitutive relations. Related efforts of the zigzag-function enrichment using higher-order 
polynomials are also highlighted. In Section 2, the principle of virtual work is presented, which 
serves as the foundation for developing finite element approximations. Also highlighted in this 
section are (i) a simple and effective method of computing highly accurate interlaminar stresses, 
(ii) the unique modeling of homogeneous plates using the full kinematics of RZT, and (iii) a 
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higher-order theory that includes odd and even zigzag kinematic terms in the inplane expansions 
and which also accounts for the thickness-stretch deformations. In Section 3, recent efforts to 
develop RZT-based finite element for laminated composite and sandwich beam, plate, and shell 
structures are discussed, focusing on their kinematic approximations, element nodal 
configurations, and modeling capabilities. In Section 4, finite element results for a laminated 
cylindrical shell undergoing large displacements are presented to demonstrate the latest RZT- 
based shell modeling capability. Finally, conclusions are presented which highlight the salient 
features of RZT and the finite elements developed on its foundation. 

2 FOUNDATION OF REFINED ZIGZAG THEORY 

Consider a multilayer composite plate of uniform thickness 2 h that is composed of perfectly bonded 
orthotropic layers, labelled with superscript (k), where k=\ ,..., A; furthermore, the plate undergoes small- 
strain deformations under static loading while exhibiting negligible inertial effects (refer to Figure 1 .) The 
inplane coordinates of the plate are defined by the vector x = (Xj , x 0 ) e S m , where S m represents the set 
of points given by the intersection of the plate with the plane z — 0 (the midplane); the symbolism 
Z e [-h, h ] denotes the domain of the through-the -thickness coordinate. 
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Figure 1. (a) Plate of uniform thickness, 2 h, subjected to transverse pressure loading, d , and edge 
tractions, { T a (a = 1,2), T_ J; (b) through-thickness notation of N -layer laminate. 


The displacement vector of any material point P(x l ,x 2 ,z) in the k th layer is defined by the three 
orthogonal components (u[ k) ,u 2 ] ,u_) that are expressed as [4-5], 

u[ k) (x, Z ) = m ( x ) + Z 9\ (x) + <f>l k) {z) ^i(x) 
u 2 \x,z) = v(x) + z 9 2 {x) + (j) 2 \z) ^ 2 ( x ) 

U z (x,z) = w(x) . 
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where u(x) and v(x) are the inplane displacements in the x 1 and x 2 coordinate directions, 
respectively, and w(x) is the transverse deflection. The symbols i9 t ( x) and 0 2 (x) represent overall 
bending rotations of a transverse material line element, that is initially straight and normal to the 
midplane, about the positive x 2 and the negative x, directions, respectively. These overall rotations 

represent linearized weighted averages of the actual nonlinear through-the -thickness displacement 
variations that are produced when the displacement variations within each layer are pronounced. 
Collectively, these three displacements and two rotations form the basis of FSDT and, as such, are viewed 
herein as the quantities that define the overall, first-approximation, coarse displacement characteristics of 
a plate. 

Refinements to the coarse displacement characteristics, and hence improved fidelity, are obtained by 
introducing additional displacement functions that model lamina-scale, or even sub -lamina-scale, 
responses adequately. In the RZT, the refined contributions to the “coarse” inplane displacements are due 

to the zigzag terms </> ( a k \z) y/ a (x) (a = 1,2) that appear in Eqs. (1). These refinement functions 
generally produce piecewise-continuous, nonlinear through-the -thickness distributions that are defined by 
the zigzag functions and their corresponding rotation (or amplitude) functions y/ a (x ) . The 

simplest form for ^>„ k \z) is given by a set of piecewise linear, continuous functions with derivatives 

(f) {k \ that are discontinuous at the layer interfaces. Herein, a comma followed by the subscript z denotes 
differentiation with respect to the through-the -thickness coordinate z. This particular class of functions is 
denoted herein by C <0) . 

The localized rotations of the line elements within each layer, which prior to deformation are straight 
and perpendicular to the reference midplane, are obtained from Eqs. (1) as 

u al ( X ’ Z ) = 9 a ( X ) + fta ¥ a ( X ) (« = 2 ) ( 2 ) 

where j3^ k) = <p k \ are constant valued within the k th layer and generally have different values for layers 
with different material properties; thus, they are discontinuous at the layer interfaces. Therefore, the 
rotations u (k \{x, z) are also uniform within the k th layer and generally vary from layer to layer. 

A key step within RZT is to identify its relationship with FSDT. This task is accomplished by 
representing the coarse rotations, 6 a (x) , as the weighted-average rotations, i.e., 

= ( a = 1 ’ 2 > 

2h (3) 

Substituting Eq. (2) into (3), results in 

f * Pa ) dz = </> ( a N) (h) - ^ (-K) = 0 (a = \,2) 

J-h (4) 

thus guaranteeing equal values of the zigzag functions on the bounding surfaces, i.e., 
= (p^P-h) (a = 1,2). To correlate the displacement of FSDT and RZT at the bounding 

surfaces, z — ±/t , it becomes immediately apparent that the top and bottom values of the zigzag functions 
must vanish identically, i.e., 
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ff\h) = $\-h) = 0 (a = 1, 2) 


(5) 


Equations (3) also imply that the following weighted-average transverse shear strains are those which 
correspond to FSDT, i.e., 

Ya = ^7 \ h . Yaz d Z = W a + (« = 2 ) 

2 h J ~ h ( 6 ) 

where henceforth the notation (•) a = d(») / dx a denotes partial differentiation with respect to the 
midplane coordinate, x a . 

It can now be readily ascertained that within RZT only the coarse kinematic variables define the 
values of the inplane displacements at top and bottom surfaces, i.e., 

u[ l) (x, - h ) = u(x) - h 9 X (x), u[ N) (x, h) = u(x) + h 6 X (x) 

u 2 \x, -h) = v(x) - h d 2 (x), u 2 r \x, h) = v(x) + h 0 2 {x) ^ 

With this insight, it becomes clear that the zigzag kinematics can indeed be regarded as local ply-level 
perturbations to the overall, coarse displacement fields. From Eqs. (7), the physical interpretations of the 

coarse variables { u , V , 0 X , 0-, j in terms of the top-surface inplane displacements u at (x) = //^ V) (x, h) 
and the bottom-surface inplane displacements u ah (x) =u^\x,—h) also become apparent, i.e., 


, , 1 r 

m(x) = 2 l 


U ]t (x) + u U) (x) | , 0 l (x) = - 1 7 r u lt (x) - u ]h (x) 

2 n L 


v(x) = - r. u 2l (x) + u 2b (x)l , e 2 (x) = - 1 - 1" u 2t (x) - u 2h (x) 
2 L J 2h L 


These same relations hold for these kinematic variables used in FSDT. 


(8) 

(9) 


For certain classes of problems, modelling advantages can be exploited by using a different set of 
primary unknown functions instead of the seven appearing in Eqs. (1). For example, in finite element 
analyses that use the sub-laminate concept [6-9], the inplane displacements of the top and bottom 
bounding surfaces of a sub-laminate are particularly useful as primary unknowns. Thus, by treating a 
laminate as a collection of contiguous sub-laminates, the RZT equations presented herein can be modified 
as follows. From Eqs. (8) and (9), it is seen that an alternate description of the RZT displacement fields is 
obtained by expressing the four coarse variables { u , V , 0^,0-, ] in terms of the four inplane 

displacements of the top and bottom surfaces { u at (x ) , u ab (x) } (a = 1,2), and then substituting the 
result into Eqs. (1). In this way, RZT is redefined in terms of the alternative set of seven displacement 
variables: { u at ( x) , u ab (x) , \// a (x) , w(x ) } (a — 1,2) which can be applied on a laminate and a sub- 
laminate level. Because this transformation of the primary unknowns is well defined, this version of RZT 
is expected to have the same analytic accuracy as the original RZT. 

Using the linear strain-displacement relations of elasticity theory and the displacement assumptions 
given by Eqs. (1), the RZT strains are given as 
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(10) 


£ ii = u ,i + + 4i Vu 

^22 ' = V ,2 + ~^2,2 + 4*2 ) Y 2,2 

Yl 2 ' = M ,2 + V ,1 + ~(^ 1,2 + ^ 2 , 1 ) + ^ W \,2 + 4 l '' W 2,1 

Yal=Ya+PaVa (« = U) 


where it is noted that the inplane strains, { fjf * , , /p ' }, are piecewise linear through the laminate 

thickness. The transverse shear strains, , are piecewise constant , i.e., they are constant within each 

material layer but discontinuous along the layer interfaces, in contrast to those of FSDT which are 
constant across the total laminate thickness. 

Each layer of the plates considered herein is presumed to be specially orthotropic with respect to a set 
of principal material coordinate directions. In addition, each set of layer principal material axes have two 
axes that reside in the plane of the layer and a third that is perpendicular to that plane. This third axis is 
parallel to the z-axis for the plate coordinates defined herein. The two principal axes that reside in the 
plane of the layer are generally noncoincident with the inplane coordinate axes of the plate. In terms of 

the plate coordinate system, the generalized Hooke's relations for the k th layer are given herein by 
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for the inplane stresses and by 
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or 


T W _ qW y(k) 


( 11 . 1 ) 


( 11 . 2 ) 


for the transverse shear stresses. The symbols and } denote the transformed elastic-stiffness 

coefficients referred to the plate coordinate system. In RZT, these coefficients are also based on the 
presumption that the transverse-normal stresses are negligible. 

The next step in formulating RZT is to provide a mathematical description of the zigzag functions. 

Thus, for the k' h material layer, which is located in the range \z, (k | )5 Z. (k] ] , the zigzag functions are given 
as (see Figure 1(b) depicting the lamination notation), 

(a = 12) (12.1) 

where 

g(k) = 


(z-z^j/h^-l e[-l,l] (* = 1,...,#) 


(12.2) 
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with the first layer, k ~ 1 , beginning at z (()) = —h , the last N layer, k — N ^ ending at z. (N) = h , and 

the * layer ending at z (/f > =z (/t _ 1) +2/7 , where z/7 denotes the ^ layer thickness. Evaluating 

Eqs. (12) at the bottom (<^ ( k) = — 1) and top = 1) surfaces of the k layer, gives rise to the 

definitions of the interface displacements 


4(,-i)=C ) (C a ' ) =-D, <t> a(k) =€\C (k) =r> {k = 


(13) 


whereas according to Eqs. (5) the zigzag functions vanish at the bottom (k = 1, <^ (1) = — 1) and top 
C k = N, C, {N) = 1) plate surfaces, i.e., (j) a{0) = (j> a(N) = 0 . 
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k= 2 

h 

l 

f 



2 r 

t ) 
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(a) Layer notation 




(b) Zigzag functions (j)\ k) (z) and (j)~, k) (z) . 


Figure 2. Notation for a three -layer laminate and (f) ik) {z.) zigzag functions defined in terms of interface 
displacements, (f) (/ik) . 
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The remaining interface displacements can be obtained from the simple expression 


0a(k) ~ 



(14) 


In [5], two alternative expressions have been derived for /3 { a k) ( k =1,..., TV — 1 ); they are given as 


(a) 


A k) - 


i 


Ak) 


a 


(k) 


h 


" h (k) 

2-i Mk) 
k= 1 xiaa 


-1. (b) /?<*> = 


-j\h (k) S<2 

ht aa 


-1 (a =1,2) 


(15) 


where Q ( Jj and (a — 1,2) denote, respectively, the diagonal transverse-shear stiffness and 

compliance coefficients of the k th layer, respectively. Note that these expressions are only slightly 
different from each other for the general lay-up case, in which the transverse-shear constitutive matrix has 

a non-zero off-diagonal term (Q^ ^0). The two definitions, however, are identical for the special, 

decoupled case for which Q [ k ' = 0 . From the predictive perspective, either form of /3 ^ ) yields excellent 

results, with the second form exhibiting slightly improved results, particularly for sandwich plates. 

Further accuracy enhancements within the present RZT methodology can be achieved by using higher- 

degree polynomial expansions for the zigzag functions </>^ 1 , that are piecewise continuous nonlinear 
functions, while keeping the basic kinematic assumptions, Eqs. (1), unchanged (e.g., refer to [10-13]). 


2.1 Variational framework 


The principle of virtual work, for the case of negligible body forces and zero shear tractions on the top 
and bottom bounding plate surfaces, may be written as 

h 


JJ( a (k)T Se (k) + r (k)T Sy (k) ) dz clS - J q SwdS 

S m — h S m 

h 


(16) 


- J J [T^ +T 2 du[ k) + T z 8u, 


dsdz = 0 


C a ~h 


where 8 is the variational operator; <7 is the applied transverse pressure attributed to the middle reference 
surface, S m ; { T a (a = 1, 2) , T } are the inplane and transverse shear tractions that are prescribed along 
the cylindrical edge surface of the laminate, S a = C a x s . The bold-face vectors appearing in equation 
(16) are defined by Eqs. (11). 

Integrating Eq. (16) across the laminate thickness, while accounting for the relationships in Eqs. (1) 
and (10), yields the corresponding two-dimensional statement of the principle of virtual work 



J s [N> m+ M ^e fc+ Q> s -^w]J5-J c F ? Au ds = 0 


(17.1) 


where the membrane stress resultants and conjugate strain measures are defined as 


={Ni,N 2 ,N 12 } = \ h _ h {a‘[\\ < T^ ,r^}dz, ^ T m = { u^v a ,u 2 +v d ) 


(17.2) 


Likewise, the bending stress resultants and conjugate strain measures are defined as 


Ml = [M v M ( , M 2 , M+, M n , M.( 2 , M{ x , 

= \ h h { Z(T U' $ k) za 22, 4 k) °22’ ZT$ ,(j{ k) T$\ 4>2 ] T$}dz 

e b = {^ 1 , 1 ’ ^ 1,1 > ^ 2 , 2 ’ V 2,2’ @1,2 + @2.1 ’ ^ 1,2 > ^ 2,1 } 


(17.3) 


The transverse shear stress resultants and conjugate strain measures are defined as 


Q! - {&, Qt e„ Qt} = J_‘ {r<» PP rg>, r«>, A<‘» r«>}& 

= {w 2 + 0 2 , y/ 2 , Wj + <9 I? 


(17.4) 


The force and moment resultants due to the prescribed tractions and their conjugate displacements are 
defined as 


F r 


u 


T 


{u,V,W,0 l ,0 2 ,l// 1 ,l// 2 } 


(17.5) 


All elements of Eqs. (17) with the superscript (j) are associated with the zigzag functions. 

The stress resultants are readily obtained in terms of the two-dimensional plate strain measures by 

integrating the expressions in Eqs. (17) through the laminate thickness, resulting in the following 
constitutive relations of RZT : 
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( 18 ) 


N m 


A B 0 



.= 

B r D 0 


Q s 


0 0 G 

. e , s . 


where A = [^y] 3x3 denotes the membrane stiffness matrix, B = [S ;/ ] 3x7 is the membrane-to-bending 
coupling stiffness matrix, D = [D (/ ] 7x7 is the bending stiffness matrix, and G = [G (/ ] 4x4 is the transverse 
shear stiffness matrix. The explicit forms of these stiffness coefficients are found in [4,5]. Note that for a 
general composite lay-up, the membrane-to-bending coupling stiffnesses are nonzero; that is, B ~ * 0 . 


2.1.1 Special case 

When the zigzag functions in the x 1 and x 7 directions are identical, 

(z) = (f> 2 k> (z) = (/) {k \z) (k = . As a result, the two twisting moments associated with the 

zigzag kinematics are given by the same expression (refer to Eqs. (17.2) and [14]) as follows 

Mt 2 =Mi x =\ h _/ k \^dz (19) 

This simplification gives rise to the following reduced form of Eqs. (17.2) 

M l ={M { M 2 ,M n ,Mf 



e b ~ {^ 1,1 ’^ 2,2 ’^ 1,2 + ^ 2,1 iWl ,! ’¥\,1 + V / 2 , l } 

and the resulting simplifications in Eqs. (18). These constitutive relations provide a clear physical 
interpretation of the new quantities associated with the zigzag kinematics. Specifically, 

} represent the bending and twisting moments due to the zigzag related cross-sectional 

distortions of the layers, whereas {y/^ l5 y/ 2 2 , (//, 2 + (/A | } are their conjugate curvatures due to the normal 

zigzag rotations, y /\ an ^ Wi ■ 


2.2 Interlaminar stresses 

An important attribute of any plate theory that is formulated to be robust at the local level is the ability to 
predict adequately interlaminar stresses that can be used to assess laminate failures. To provide adequate 
estimates of interlaminar stresses, it is often customary to integrate the three-dimensional equilibrium 
equations of elasticity theory. In this approach, the interlaminar transverse shear stresses are obtained by 
integrating the x 1 and x 2 derivatives of the inplane stresses 
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( 21 . 1 ) 



v°l 1,1 

(CT {k) 

1 ° 22,2 



dz 


whereas the interlaminar transverse normal stresses are obtained by integrating the derivatives of the 
transverse shear stresses as 


Z 



-h 



(21.2) 


Although this post -processing procedure has been widely advocated for use with FSDT and many higher- 
order theories, the adequacy of this scheme can be seen to be directly linked to (a) the accuracy of the 
inplane stresses obtained from constitutive relations of the underlying theory, and (b) the accuracy of the 
derivatives of these stresses (or second derivatives of the kinematic variables, refer to Eq. (10)), with this 
issue particularly relevant when finite element analyses are performed. 

Regarding issue (a), it has been observed that FSDT and the vast majority of available structural 
theories often underestimate the inplane stresses; this is especially the case when highly heterogeneous 
and relatively thick laminates are analyzed, including the sandwich construction. Relative to issue (b), 
reasonably good improvements in estimating derivatives of the inplane stresses can be achieved by using 
smoothing techniques, e.g., [15]. 

Application of Eqs. (21), using RZT inplane stresses, has been shown to be highly effective, yielding 
accurate interlaminar stresses that are comparable to those predicted by three-dimensional elasticity. This 
robustness exists because the piecewise linear inplane stresses of RZT are consistently accurate, even for 
highly heterogeneous and relatively thick composite and sandwich laminates [1-5]. 


2.3 Modeling of homogeneous plates 

The key property of a zigzag function is that it vanishes identically when the material has homogenous 
transverse shear properties across the total thickness. Consequently, when the zigzag terms vanish 
identically, the kinematic assumptions given by Eqs. (1) revert to those of FSDT, which require transverse 
shear correction factors. Tessler et al. [5] have found a simple and effective way to use the full power of 
RZT’s kinematic field to model homogenous plate (and beam) problems, without the need for transverse 
shear correction factors. In this approach, a homogeneous plate is modelled as a multilayer heterogeneous 

laminate, in which the transverse shear moduli, G^(a = 1,2) , vary only slightly from layer to layer, i.e., 


G^=G a3 (l + e (k) ) (cc = 1,2; k = (22.1) 

where (i a = 1, 2) denote the constant values of the shear moduli, and where 

S' j (^(fc-i) + Z(fc) + Z(k-l)Z(k) ) ^ 1 ( s ^ 1) (22.2) 

is a layerwise, dimensionless coefficient that is a function of the position of the & th layer within the 
laminate thickness. Thus, since S (k> <SC 1 (e.g., by setting .v = 10 4 ), the values of are only slightly 
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perturbed compared to the constant shear-modulus values of the corresponding homogeneous plate. 

Thus, in this approach, the homogeneous plate is replaced with a corresponding plate with an 
infinitesimal degree of heterogeneity, which for all practical purposes is homogeneous. This approach is 
effective even when the deviation from the constant value of the shear modulus of the homogeneous plate 
is less than 1/100 of a percent. The results in [5] have shown that as the number of material layers 
increases within this infinitesimally heterogeneous laminate, the solution approaches that of a 
homogeneous plate, exhibiting the well-known parabolic distribution of the transverse shear stresses 
across the plate thickness, as well as the cubic distribution of the inplane stresses, with the latter being 
particularly evident in relatively thick plates. Remarkably, the solution itself finds the correct transverse 
shear stress distributions, albeit in a piecewise constant manner, thus validating the lack of necessity for 
the transverse shear correction factors that are commonly used within FSDT and even within some 
theories of higher order. 


2.4 Inclusion of transverse normal deformations 

The RZT presented herein incorporates FSDT as the underlying baseline plate theory that models the 
overall, coarse plate behavior. Thus, it follows that higher order RZT theories can be formulated by using 
other underlying baseline plate theories with the zigzag kinematics used herein or with zigzag kinematics 
of higher fidelity. For example, Barut et al. [14] combined an underlying { l,2}-order baseline theory, that 
uses an overall parabolic transverse displacement assumption, with piecewise quadratic RZT-based zigzag 
inplane displacements to obtain a higher-order refined zigzag plate with eleven kinematic variables; four 
more than the RZT presented herein. The displacement assumptions for this theory, labeled as RZT 12 ' 21 , 
are given as 


u[ k \x, Z) = u(x) + Z%(x) + tf k) (z)yr n (x) + | # } (z)^ 12 (x) 

h 


u ( 2 k \x,z) = v(x) + z0 2 ( x ) + </> 2 ) (z)i/f 21 (x) + j<f>y\z)ii/ 22 (x) 

h 


(23) 


u (X,z) = w(x) + f Wj(x) + 
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v /r 3y 


w 2 (x) 


where <f>^ are the piecewise linear zigzag functions within the A:' 1 ' layer, adopted from the original RZT 
formulation presented herein; y/ a] (x) and y/ al {x) (a = 1,2) are the four zigzag amplitudes that permit 
both non-symmetric, i// a] (x), and symmetric, x// a 2 (x), inplane zigzag deformation modes. In addition, 
an average transverse normal stress, cr„ , is independently assumed as a cubic function through the 
laminate thickness, in the form proposed by Tessler [16]; that is 


o- 77 (x,z) = o\ 0 (x) + cr 7l (x) 


'zl' 
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This representation of the normal through-the-thickness stress yields a parabolic expression for its 
derivative <x_ Z; _ (x, z) , given by 


12 



( 


(24.2) 


<7zz,z( X ’Z) 


^-l( x ) 

h 


±h) = 0 


.2 ^ 



(24.3) 


Equations (24.3) constitute the exact conditions on a.. in the sense of three-dimensional equilibrium 
equations of elasticity theory for the case of zero-valued transverse-shear tractions on the bounding 
surfaces. The cr_ 0 (x) and cr,j (x) terms in Eq. (24) are expressed in terms of kinematic variables of the 

theory by way of the least-square transverse strain compatibility relations originally proposed in [16]. 

Thus, in this new theory, both transverse shear and transverse normal (thickness-stretch) deformations 
are included. The present theory models accurately the in-plane and transverse stress components through 
the thickness. The addition of the even (symmetric) inplane zigzag modes, which are not included in the 
seven-variable RZT, contributes to some improvements of the inplane displacement, strain and stress 
response, with further improvements in the interlaminar transverse shear stresses. Furthermore, RZT {2 ' 21 
appears to be an excellent candidate for developing efficient and accurate C°-continuous finite elements. 


3 FINITE ELEMENT APPROXIMATIONS 

The variational statement given by Eq. (17.1) can be integrated by parts to produce a consistent set of 
equilibrium equations and boundary conditions, resulting in seven partial differential equilibrium 
equations in terms of seven kinematic variables [4-5]. The equilibrium equations can be solved exactly or 
approximately depending on the complexity of the material lay-up, boundary conditions, and loading. 
Refer to [1-5] for details of analytic and approximate solutions of simply supported and cantilevered 
beams and plates made of laminated composite and sandwich construction. 

Alternatively, to enable large-scale analysis of complex aerospace structures, the variational statement, 
Eq. (17.1), can be discretized using finite elements. As in FSDT, the strain measures within RZT are 
given as first-order partial derivatives. The direct implication is that computationally efficient C°- 
continuous beam, plate, and shell finite elements can be developed, and that the legacy finite element 
technology that has been developed for FSDT can also be used with RZT. 

3.1 RZT beam elements 

Recently, Gherlone et al. [17,18] and Onate et al. [19] derived planar beam finite elements based on the 
RZT beam formulation in Tessler et al. [1]. When the beam deformations are restricted only to the Xj - z 
plane, the second equation in Eqs. (1) is identically zero, and the theory reduces to the two displacement 
components, u\ and w, . Consequently, there remain four independent kinematic variables, 
u(x l ), vv(X| ), O, lx, ), and i// l (x, ) . that describe the membrane, bending, transverse shear, and zigzag 
deformations. 

Gherlone et al. [17,18] derived several low-order RZT beam finite elements with the aim of achieving 
the best compromise between accuracy and computational efficiency. The four kinematic variables use 
the anisoparametric (aka interdependent ) interpolations, where the polynomial degree of vv(X| ) is one 

order higher than those approximating the u(x ] ) , 0, (x, ) , and I//, (x, ) variables. Such interpolation 
strategy enables free of shear locking element behaviour for slender beams. With an initial assumption of 
a parabolic distribution for wix ] ) , the authors explored several shear-related constraint conditions that 

gave rise to two-node elements and a coupled form for the vv(Xj) interpolation. The constraint condition 
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requiring a constant variation of the transverse shear force gave rise to a remarkably accurate two-node 
beam element. For further details of this formulation and for elastostatic solutions of simply supported 
and cantilevered beams of various slenderness ratios and lamination properties, the reader is referred to 
[18]. 

Onate et al. [19] derived a two-node, RZT beam finite element using linear, isoparametric 
interpolations for the four kinematic variables, u(x ] ) , w(Xj ) , 9 { (x j), and y/ } (x l ) . To achieve an 

element capable of modeling slender beams without incurring the shear locking effect, the authors used 
reduced integration (a one point Gaussian quadrature rule) of the transverse-shear strain energy of the 
beam element. One of the most interesting numerical solutions presented was for a laminated beam with 
an imbedded delamination that was modeled as a compliant layer. It was shown that the RZT beam 
element is computationally efficient and is capable of high-fidelity modeling of imbedded delaminations. 

3.2 RZT plate and shell elements 

Versino et al. [20,21] developed six- and three-node triangular RZT-based plate finite elements. Adopting 
linear shape functions for the in-plane displacements, bending rotations, and zigzag amplitudes, and a 
quadratic shape function for the transverse deflection, (i.e., using the Tessler-Hughes anisoparametric 
interpolation strategy [22]), the element interpolations in terms of the linear area-parametric coordinates 
L ( have the form 


3 3 6 

u(x v x 2 ) = Y J u i L i, v(x 1 ,x 2 ) = Y J v i L i , w(x 1 ,x 2 ) = Y J w k P k 

i= 1 i= 1 k= 1 

3 3 

d\ (*1 , *2 ) = Z i L i ’ (x l ’ *2 ) = Z 9 2i L i (25) 

i=l i=l 

3 3 

V\ Ul > *2 ) = Z V'lA ’ ^2 (*l ’ X i) = Tj ^2 i L i 

i'=l i=l 

where i e {1, 2, 3} is an index ranging over the three comer nodes; k e {1, m, 2 , 2, m 2i , 3, m 31 } ranges over 

the corner and mid-edge nodes; and P k defines the standard set of quadratic shape functions (see the 

nodal configurations in Figure 3.) The element topology includes six nodes; however, the mid-edge 
nodes have only degrees-of-freedom (dof) associated with the transverse deflection. Using these 
interpolations, the six-node element, called unconstrained, is readily derived by introducing Eqs. (25) in 
the variational principle, Eq. (17), while carrying out the necessary matrix and variational operations to 
obtain the element stiffness equations. 

In addition, the authors derived a three-node, constrained anisoparametric element whose interpolation 
functions for the deflection variable in Eqs. (25) are constrained by suitable edge constraints, analogous to 
those which have been explored by Tessler and Hughes [22] for plates modeled with FSDT, and by 
Gherlone et al. [17,18] for beams modeled with RZT. The explicit edge-constraint procedure enforces the 
mid-edge w dof to be dependent on the corner-node dof of the rotation variables of the corresponding 
edges. By imposing one-dimensional constraint relations, requiring the two shear strain measures to be 
constant along the element edges, there results a coupled-form, parabolic deflection that is compatible 
along the element edges given by [21] 

3 3 

w(Xi ,x 2 ) = X L ;W,+X[K + cx Pn ) K + ( ^ + Wn ) U, ] (26. 1 ) 

1=1 1=1 
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with 


L U~^( b k L j h , L k)’ ^2, - 2 ( a j L k a k L j ) (26 2) 

a i= X lk~ X lj’ b i= X 2j~ X 2k 

where the subscripts are given by the cyclic permutation of i e { 1, 2, 3} , / e { 2, 3, 1} , and /< e {3, 1, 2} ; 
and where c is either 0 or -1, depending on the constraint strategy used in the element formulation. 



Figure 3. Six-node (unconstrained) and three -node (constrained) RZT-based anisoparametric plate 

elements. 

An extensive numerical study was carried out in [21] on various laminated composite and sandwich 
plates undergoing elasto-static deformations. The results of this study demonstrate that the RZT-based 
anisoparametric elements possess superior accuracy and excellent convergence characteristics over a wide 
range of the span-to-thickness ratio. In addition, it was shown that the three -node elements provide the 
best compromise between computational efficiency and accuracy. Furthermore, the RZT models provide 
superior through-the -thickness predictions of the inplane and interlaminar stresses over comparable FSDT 
models, especially for relatively thick and highly heterogeneous laminated plate and sandwich plate 
constructions. 

The latest finite element implementation of RZT, by Versino [23], is focused on developing 
robust anisoparametric shell finite elements that include the drilling dof (i.e., the dof associated 
with a rotation about the normal to the element’s planar surface) and geometrically nonlinear 
deformations. In the next section, the application of RZT to an elastic cylindrical shell 
undergoing large displacements under quasi-static loading is demonstrated. 


4 NONLINEAR DEFORMATIONS OF A LAMINATED CYLINDRICAL SHELL 

In this demonstration problem, a cylindrical shell is subjected to a quasi-static concentrated force, F, 
applied at point C (Figure 4). The global shell dimensions are L=254 mm, r=2540 mm, and 9=0.1 rad. 
The shell wall is a three-layer laminate with a total thickness 2/7=24 mm and with each layer having the 
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same thickness. The prescribed boundary conditions are such that the straight edges are fully clamped and 
the curved edges are free. The mechanical material properties of isotropic material layers A and B are 
summarized in Table 1, and the laminate stacking sequence is given by [A/B/A], Note that the elastic 
modulus of the middle layer is two orders of magnitude less than those of the top and bottom layers. This 
aspect, typical of sandwich construction, presents a significant challenge for any structural theory, since 
the through-the-thickness distributions of the inplane displacements exhibit interface slope discontinuities 
that are not accounted for in most structural theories. 

To establish a viable reference solution for this problem, a high-fidelity, three-dimensional 
geometrically nonlinear FEM analysis was performed using the ABAQUS commercial code [24], 
Depicted in Figure 4(a) is an FEM model which is arrived at from a convergence study. The model is 
based on C3D20 brick elements, discretizing a symmetric quadrant of the shell, using a 64x64x6 mesh 
(64 elements along each edge and 2 elements across the thickness of each layer). In addition, two shell- 
based geometrically nonlinear solutions were obtained: (a) an ABAQUS solution using the S3 (three- 
node, FSDT) shell elements, and (b) an RZT-based solution using RZT3 - a three-node anisoparametric 
element with drilling dof. Both shell models are based on the 32x32 mesh subdivisions that span the 
shell’s symmetric quadrant, as shown in Figure 4(b). 

Comparisons of the three nonlinear FEM solutions are shown in terms of the load-deflection curve at 
point C in Figure 5. It is evident from these results that the three-dimensional FEM solution and the 
corresponding RZT3 solution are in very close agreement over the entire range of the applied loading. In 
contrast, the S3-model predictions are considerable less accurate, particularly over the range of larger 
displacements. 


Table 1. Mechanical material properties for isotropic layers. 


Material 

Young’s 
modulus, 
E (MPa) 

Poisson 

ratio, 

V 

A 

31 x 10 2 

0.3 

B 

31 
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Figure 4. Finite element discretization for cylindrical shell. 



Figure 5. Applied force vs. deflection at point C. 


5 CONCLUSIONS 

In this paper the latest advances in the development of a shear-deformable, Refined Zigzag Theory (RZT) 
for homogeneous, laminated composite, and sandwich beams and plates have been reviewed and the 
salient features highlighted. This theory maintains four kinematic unknowns for beams and seven for 
plates, regardless of the number of layers, material composition, or layer stacking sequence. In addition, it 
has been shown that the theory can be expressed in terms of displacement quantities that are amenable to 
sub-laminate modelling techniques. A higher-order plate theory that incorporates the RZT presented 
herein has also been described, which includes the thickness-stretch deformations in addition to the 
transverse shear deformations, and has eleven kinematic variables. This higher -order theory illustrates 
how the basic RZT presented herein can be adapted to obtain results with a very high degree of fidelity. 

Analytic and computational aspects of the RZT and several new finite elements for beams, plates, and 
shells have also been highlighted. The unique characteristics of RZT and its finite elements is that they 
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permit efficient modelling of a wide range of problems including homogenous beams and plates, 
composite laminates, sandwich construction, and compliant-layer delamination modelling. Such 
solutions do not require shear correction factors or any increase in model complexity or computational 
effort. In addition, RZT-based finite elements are amenable to much of the legacy finite element 
technology, especially that for FSDT-based elements. Because of these attributes, RZT has considerable 
promise for becoming the theory of choice for many practical applications including the computationally 
challenging problems of progressive damage modelling in composite structures. 
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